#   LibAut
#       Episode criteria sensitivity 
#       Appendix B, Table B1, Figure B1 and B2
#   Juraj Medzihorsky
#   2021-08-11

library(parallel)
options(mc.cores=3, stringsAsFactors=F)

gr <- (1+sqrt(5))/2

#   sessionInfo()
##  R version 3.6.3 (2020-02-29)
##  Platform: x86_64-pc-linux-gnu (64-bit)
##  Running under: Ubuntu 16.04.7 LTS
##
##  Matrix products: default
##  BLAS/LAPACK: /opt/OpenBLAS/lib/libopenblas_haswellp-r0.3.9.so
##
##  locale:
##   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
##   [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_GB.UTF-8       LC_NAME=C
##   [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
##  attached base packages:
##  [1] parallel  stats     graphics  grDevices utils     datasets  methods   base
##
##  other attached packages:
##  [1] nvimcom_0.9-28 colorout_1.2-0
##
##  loaded via a namespace (and not attached):
##  [1] compiler_3.6.3 tools_3.6.3

#------------------------------------------------------------------------------
#   data loading

code_dir <- getwd()
data_dir <- gsub('code$', 'data', code_dir)
figs_dir <- gsub('code$', 'figures', code_dir)

setwd(data_dir)
v <- readRDS('V-Dem-CY-Full+Others-v10.rds')
e <- ert <- read.csv('eps.csv')
along <- readRDS('episode_sensitivity_20210810.rds')
setwd(figs_dir)

ert$is_libaut <- with(ert, (dem_ep_prch%in%1) & (dem_ep_ptr%in%0))
v <- subset(v, year>=1899)

#   this is in Table B.1
val_grid <-
    expand.grid(
        si=seq(0.005, 0.015, 0.001),
        ci=seq(0.05, 0.15, 0.01),
        yt=seq(0.01, 0.05, 0.005),
        ct=seq(0.05, 0.15, 0.01),
        to=seq(3, 7, 1)
    )

#------------------------------------------------------------------------------
#   data transforms 

along$prop <- along$Freq/nrow(val_grid)
atab <- xtabs(prop ~ uid + dem_ep_outcome, data=along)
class(atab) <- 'matrix'
adf <- data.frame(atab)

a <- data.frame(matrix(nrow=nrow(atab), ncol=0))
a$country_text_id <- sapply(rownames(atab), function(i) strsplit(i, '_')[[1]][1])
a$year <- sapply(rownames(atab), function(i) as.numeric(strsplit(i, '_')[[1]][2]))

a$p_suc <- adf[, 'X1']
a$p_fai <- rowSums(adf[,c('X2','X3','X4')], na.rm=T)
a$p_cen <- adf[,'X6']
a$p_epi <- rowSums(adf, na.rm=T)
a$p_out <- 1-a$p_epi

id_vars <- c('country_text_id', 'year')

vs <- merge(v[,id_vars], a, by=id_vars, all.x=T)
vs$p_suc[is.na(vs$p_suc)] <- 0
vs$p_fai[is.na(vs$p_fai)] <- 0
vs$p_cen[is.na(vs$p_cen)] <- 0
vs$p_epi[is.na(vs$p_epi)] <- 0
vs$p_out[is.na(vs$p_out)] <- 0

dat <- merge(vs, ert)


#   clumsy color lightener
lighten <-
    function(color, bw=0.5)
    {
        kol <- as.numeric(col2rgb(color))
        num <- pmax(0, pmin(255, bw*rep(255,255,255) + (1-bw)*kol))
        nkl <- rgb(num[1], num[2], num[3], max=255)
        return(nkl)
    }


#   nature color blind friendly palette again
mks <- c(rgb(  86, 180, 233, max=255),  #   sky blue
         rgb( 230, 159,   0, max=255),  #   orange
         rgb(   0, 158, 115, max=255),  #   bluish green
         grey(0.5))  #   grey
mka <- sapply(mks, lighten)


#------------------------------------------------------------------------------

pdf('Figure_B1.pdf', 5, 6)
#
par(mfrow=c(4,1), mar=c(0.5,0.0,2,0.0), oma=c(1,0,2,0))
a1 <- seq(0, 1, by=0.2)
#
tf <- dat$p_epi[(dat$dem_ep_outcome%in%c(1,5))&(dat$is_libaut)]
dn <- density(tf, from=0, to=1)
plot(dn, lwd=gr^1, col='white', frame=F, yaxt='n', xlab='', ylab='', main='', xaxt='n')
polygon(c(dn$x[length(dn$x):1], dn$x), c(rep(0, length(dn$x)), dn$y),
        col=mka[1], border=mks[1])
pu <- par('usr')
text(pu[1], pu[4],
     paste0('In successful episodes\nN = ', format(length(tf), big.mark=',')),
     pos=4, col=mks[1], cex=1, xpd=T)
#
tf <- dat$p_epi[(dat$dem_ep_outcome%in%c(2,3,4))&(dat$is_libaut)]
dn <- density(tf, from=0, to=1)
plot(dn, lwd=gr^1, col='white', frame=F, yaxt='n', xlab='', ylab='', main='', xaxt='n')
polygon(c(dn$x[length(dn$x):1], dn$x), c(rep(0, length(dn$x)), dn$y),
        col=mka[2], border=mks[2])
pu <- par('usr')
text(pu[1], pu[4],
     paste0('In failed episodes\nN = ', format(length(tf), big.mark=',')),
     pos=4, col=mks[2], cex=1, xpd=T)
#
tf <- dat$p_epi[(dat$dem_ep_outcome%in%6)&(dat$is_libaut)]
dn <- density(tf, from=0, to=1)
plot(dn, lwd=gr^1, col='white', frame=F, yaxt='n', xlab='', ylab='', main='', xaxt='n')
polygon(c(dn$x[length(dn$x):1], dn$x), c(rep(0, length(dn$x)), dn$y),
        col=mka[3], border=mks[3])
pu <- par('usr')
text(pu[1], pu[4],
    paste0('In censored episodes\nN = ',format(length(tf), big.mark=',')),
    pos=4, col=mks[3], cex=1, xpd=T)
#
tf <- dat$p_epi[!(dat$is_libaut)]
dn <- density(tf, from=0, to=1)
plot(dn, lwd=gr^1, col='white', frame=F, yaxt='n', xlab='', ylab='', main='', xaxt='n',ylim=c(0,40))
axis(1, a1, rep('', length(a1)), lwd=0.38, tck=-3.82*1e-2)
axis(1, lwd=0, line=-0.618)
polygon(c(dn$x[length(dn$x):1], dn$x), c(rep(0, length(dn$x)), dn$y),
        col=mka[4], border=mks[4])
pu <- par('usr')
text(pu[1], pu[4],
     paste0('Not in liberalization episodes\nN = ',format(length(tf), big.mark=',')),
     pos=4, col=mks[4], cex=1, xpd=T)
#
mtext('Country-Year Inclusion Probabilities', outer=T, font=1, cex=1,line=0.618)
#
dev.off()


#------------------------------------------------------------------------------
#   functions to make the ternary triangle for the legend

#   ternary triangle grid
toter <-
    function(v)
    {
        x <- v[1]
        y <- v[2]
        z <- v[3]
        c1 <- (1/2)*(2*y+z)/(x+y+z)
        c2 <- (sqrt(3)/2)*z/(x+y+z)
        return(c(x=c1, y=c2))
    }

#   color mixer version 1, for the legend
tokol1 <-
    function(v)
    {
        if (sum(v)>0) {
            vs <- c( 86, 180, 233)*v[1]
            vf <- c(230, 159,   0)*v[2]
            vc <- c(  0, 158, 115)*v[3]
            vk <- vs+vf+vc
            return(rgb(vk[1],vk[2],vk[3],max=255))
        } else if (sum(v)==0) {
            return(grey(0.5))
        }
    }

#   color mixer version 2, for the country-years
tokolX <-
    function(v)
    {
        vs <- c( 86, 180, 233)*v[1]
        vf <- c(230, 159,   0)*v[2]
        vc <- c(  0, 158, 115)*v[3]
        vk <- (vs+vf+vc)*sum(v) + (1-sum(v))*c(123,123,123)
        return(rgb(vk[1],vk[2],vk[3],max=255))
    }


#   color legend content
n <- 1e1
s <- (0:n)/n
g <- expand.grid(s, s, s)
g <- g[rowSums(g)==1, ]
names(g) <- c('s', 'f', 'c')
xy <- t(apply(g, 1, toter))
k1 <- apply(g, 1, tokol1)


#------------------------------------------------------------------------------
#   color the country-years

#   p_detX is any type of episode
d <- dat
d$w_p_ep   <- dat$p_epi
d$w_p_nope <- dat$p_out
d$w_p_succ <- dat$p_suc
d$w_p_fail <- dat$p_fai
d$w_p_cens <- dat$p_cen

#   add country-year color
d$w_kol <- apply(d[, c('w_p_succ', 'w_p_fail', 'w_p_cens')], 1, tokolX)

#   summaries
dl <- split(d, f=d$country_text_id)
pp <- sapply(dl, function(i) sum(i$w_p_ep))
dl <- dl[pp>=1e-3]

#------------------------------------------------------------------------------
#   plot helpers for Figure B.2

#   time-series lines
xxline <- 
    function(i) 
    {
        if ((xx$year[i]-xx$year[i-1])==1) {
            lines(xx$year[i-1:0]-0.0, xx$v2x_polyarchy[i-1:0], col='white', lend=0, lwd=xx$w_lwd[i-1])
            lines(xx$year[i-1:0]-0.0, xx$v2x_polyarchy[i-1:0], col=xx$w_kol[i-1], lend=0, lwd=xx$w_lwd[i-1])
            if (xx$w_p_ep[i-1]==0) {
                lines(xx$year[i-1:0]-0.0, xx$v2x_polyarchy[i-1:0], col=grey(0.00), lty=3, lwd=gr^-2)
            }
        }
    }

#   background rectangles based on found episodes
xxrect <-
    function(i, yv, xx)
    {
        is_lib <- xx$is_libaut[i-1]
        if (!is_lib) {
            mykol <- 'white'
        } else if (is_lib&(xx$dem_ep_outcome[i-1]%in%c(1,5))) {
            mykol <- mka[1]
        } else if (is_lib&(xx$dem_ep_outcome[i-1]%in%c(2,3,4))) {
            mykol <- mka[2]
        } else if (is_lib&(xx$dem_ep_outcome[i-1]%in%c(6))) {
            mykol <- mka[3]
        } 
        if ((xx$year[i]-xx$year[i-1])==1) {
            rect(xx$year[i-1], yv[1], xx$year[i], yv[2], col=mykol, border=mykol)
        } 
    }


#------------------------------------------------------------------------------
#   Figure B.2

pdf('Figure_B2.pdf', 30, 22)
#
par(mfcol=c(14,13), mar=c(0,0,0,0), xpd=F)
#
for (i in 1:length(dl)) {
    xx <- dl[[i]]
    plot.new()
    plot.window(xlim=c(1900,2020), ylim=c(0,1))
    box(lwd=gr^-4, col=grey(0.8))
    invisible(sapply(2:nrow(xx), xxrect, yv=c(0,1), xx=dl[[i]])) 
    rect(1900, 0, 2020, 1, col=grey(1, 0.75), border=grey(1, 0.75))
    invisible(sapply(2:nrow(xx), xxline)) 
    mtext(names(dl)[i], side=3, line=-1.62, col=grey(0.5))
}
#
par(mar=c(0,1,0,1)*4, xpd=T)
###
plot.new()
plot.window(xlim=c(-1,1), ylim=c(0,1)+c(-1,1)*1e-1, asp=1)
    mtext('Inclusion Probability', 3, line=-1.62*2)
    ps <- rev(seq(0, 1, 2e-1))
    plwd <- gr^-2 + ps * gr^1
    plty <- rep(1, length(ps))
    legend('center', bty='n', ncol=2,
           col=rep(grey(0.5), length(ps)),
           cex=rep(1.382, length(ps)),
           lwd=plwd, lty=plty, legend=round(1e2*ps))
###
plot.new()
plot.window(xlim=c(0,1)+c(-1,1)*2e-1, ylim=c(0,1)+c(-1,1)*0e-1, asp=1)
    mtext('Episode Type', 3, line=-1.62)
    points(xy[,1], xy[,2], pch=16, col=k1, cex=1.9)
    text(0.5, max(xy[,2])+gr^-6, 'Censored', pos=3, cex=1.382)
    text(-0.0, 0-gr^-6, 'Successful', pos=1, cex=1.382)
    text(+1.0, 0-gr^-6, 'Failed', pos=1, cex=1.382)
#
dev.off()


#   SCRIPT END